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A number of simple pair interaction potentials of the carbon dioxide molecule are in- 
vestigated and found to underestimate the magnitude of the second virial coefficient 
in the temperature interval 220 K to 448 K by up to 20%. Also the third virial coef- 
ficient is underestimated by these models. A rigid, polarizable, three-site interaction 
potential reproduces the experimental second and third virial coefficients to within a 
few percent. It is based on the modified Buckingham exp-6 potential, an anisotropic 
Axilrod- Teller correction and Gaussian charge densities on the atomic sites with an 
inducible dipole at the center of mass. The electric quadrupole moment, polarizabil- 
ity and bond distances are set to equal experiment. Density of the fluid at 200 and 
800 bars pressure is reproduced to within some percent of observation over the tem- 
perature range 250 K to 310 K. The dimer structure is in passable agreement with 
electronically resolved quantum-mechanical calculations in the literature, as are those 
of the monohydrated monomer and dimer complexes using the polarizable GCPM 
water potential. Qualitative agreement with experiment is also obtained, when quan- 
tum corrections are included, for the relative stability of the trimer conformations, 
which is not the case for the pair potentials. 
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I. INTRODUCTION 



For carbon dioxide (CO2), being an important industrial chemical, numerous interaction 
potentials (IPs) have been proposed, surpassed perhaps only by water in the amount of 
computer attention it has attracted among the molecular models. There are the parametric 
fits to the ab initio potential energy surface of the dimer, such as the the IPs of Steinebrun- 
ner and coworkers,- of Bukowski and colleagues,- and of Bock et al— The most successful 
and widely used IPs, however, have been fitted against bulk properties known from obser- 
vation, such as the vapor-liquid equilibrium^ - - (VLE) or the crystal lattice parameters-, 
and still others against experimental properties of the dilute gas, such as the second virial 
coefficient.-^ One may ask why this is so, but the ab initio IPs employ a great number of 
fitting parameters, not always of clear physical origin and, even so, interfacing them with 
other IPs, for instance when studying mixtures, is technically difficult because suitable com- 
bining equations are not known for the many parameters. Moreover, this problem is not 
unique to the ab initio IPs: also some empirical IPs use truncated series expansions' 1 ^ for 
both angular and radial functions in a way which makes it difficult to interface them with 
IPs of radial site-site interactions. Hence, such IPs can find little use outside simulations of 
the neat liquid. On the other hand, great success has been enjoyed by the simple site-site 
interaction formula of one Lennard- Jones interaction center, and one atomic point charge, 
centered on every atomic site^^ 1 ^. Restricting ourselves to the rigid models, these can all be 
regarded as descendants of the original IP due to Murthy, Singer and MacDonald (MSM). 
These IPs have the very appealing property that they can be readily interfaced ("mixed") 
with existing force fields, many of which share the exact same mathematical form. 

Nature is not kind enough, however, to allow such a simplified description of the CO2 
molecule at no cost. Even though very good experimental agreement for a wide variety of 
properties is obtained by the simple, rigid model with three Lennard- Jones centers and one 
point quadrupole developed by Merker et al.,— like the EPM-2 model,— it suffers from exper- 
imental disagreement in more "basic" properties such as the carbon-oxygen bond distance. 
With respect to the VLE envelope, the most successful such simple MSM-type model to 
date is the IP due to Zhang and Duan 6 (ZD). Nevertheless, despite its high accuracy in this 
property, I report in this Note that the microstructure of the dimer, and the temperature 
dependence of the second, B 2 (T), and third, B 3 (T), virial coefficients are of much poorer 
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quality. It should be pointed out, also, that the results presented in Ref. |6| have been called 
into question.— 

One striking omission from the published CO2 IPs is that of many-body effects. That 
the CO2 molecule lacks an electric dipole moment may have dissuaded investigators from 
looking in this direction, but in the work leading up to this Research Note, extensive trials 
indicated that it is not possible with a non-polarizable IP of MSM-type to simultaneously 
fit B2(T) while keeping the experimental agreement of the VLE envelope, at least not if 
the experimental bond distance and electric quadrupole moment are to be kept intact. The 
decision was then made in favor of a polarizable IP, to be described in this Note, but the 
extra cost that the high-resolution solution of the electric field equations incurs, even for 
a single polarization site, make simulations of the vapor-liquid coexistence prohibitively 
expensive for parametrization purposes. Instead, the temperature-dependence of the fluid 
density at 200 bar was used as a test for the many-body (concentrated phase) and B2(T) for 
the two-body (dilute phase) properties. Further developments then introduced a three-body 
dispersion interaction of Axilrod- Teller type,— and B$(T) as a parametrization target. As 
a test of the validity of the IP, numerous other properties are calculated without input into 
the parametrization procedure. The model introduced in this work goes by the moniker 
of Gaussian Charge Polarizable Carbon Dioxide (GCPCDO). This is because of its great 
similarity to the highly successful GCPM water— from which it borrows most of the essential 
equations. 

This Note is organized as follows. First, in Section [XT] a description of the mathematical 
form of the IP is given, and also the details of the calculations and the targeted properties of 
the parametrization. Then in Section 1111} the results are presented and discussed. Finally, 
a brief recapitulation of the main points is given in Section IIV1 

II. MODEL AND COMPUTATIONAL DETAILS 

A. Electrostatic interaction 

We compute the interaction between partial charge q a and partial charge q@ at a separa- 
tion of r a p, through the formula 

Mq(^Q/3) = -^-^(?"a/3, T cn Tp) (1) 
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Here rj(r a p, T a ,Tp) is a function that assures that the electrostatic interactions remain finite 
at all separations; for large r a p it approaches unity. Physically, this corresponds to charges 
distributed over a finite volume in space and like this we avoid the singularity of the potential 
at zero charge separation. For point charges, rj is identical to unity and Eq. (CQ) reduces to 
the classical Coulomb law. We choose the following form for the ^-function, 



which corresponds the physical case of two Gaussian charge distributions of standard devi- 
ations r a and Tp interacting with each other.— 

For now, we shall not concern ourselves with the general case of many-body interaction, 
as given by both "static" and "dynamic" electron correlation, but exclusively take care of 
that "static" correlation from the average electric field around each molecule, i. e. electronic 
induction effects. Furthermore, we do not carry this analysis beyond the dipole induction, 
i. e. we consider only the gradient of the electric potential. An induced dipole is hence 
positioned at the center of mass of molecule I and is given by 



where Ei is the electrostatic field at that point, a± the polarizability perpendicular to 
the molecular axis and an that parallel to the same. Both a± and otii are known from 
theory^ and their average agree in magnitude with that ascertained in experiment,— 1 ^ but 
the interpolation between them has been chosen merely for convenience. 9 is the angle 
between the molecular axis and the direction of Ei, the electric field at site /. Because of the 
uncertainty in the precise form of the polarization matrix, and the approximations involved 
with the rigid rotor, the polarizabilities have been rounded to only two significant digits. In 
its turn, the electric field is computed as the sum of the contributions due to the permanent 
charges and that due to the other induced dipoles, 





jit = Ei (a± + | cos(0)|(a|| - a±)) 



(3) 




where Ef is the electric field at the center of molecule i due to the all the charges on the 
other molecules, and Ef is the electric field at the same point, but due to all the other 
induced dipoles. These are given by 
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and 

4 M = X) T wZm (6) 

where the tensor Ti m is obtained from the Hessian of equation (pQ), with - following Paricaud 
et alM - the charge width parameter of the molecular center equal to that of the center 
atom. If the dipoles are converged, the extra energy of interaction due to the polarization 
is given by,— 

N 



^X>"4 q (7) 



C/poi = 2 

i=i 

B. Dispersion interaction 

First of all, the modified Buckingham exp-6 potential is adopted between atoms a and 



0, 

6 f / r \ 1 (a 



1 - 6/laf: 



exp - LJ } - (- 



(8) 



.la/3 I V °" aPJ ) \T. 

This represents the pairwise additive part of the dispersion interaction as well as the steric 
repulsion between atoms. Also for this interaction, a hard-core is introduced at r a p = 0.57a a /3 
to avoid the spurious behavior of this potential at short range. This is the same hard-core 
cutoff used by Paricaud and coworkers^ in the GCPM water model. Investigations indicated 
that the results are not very sensitive to the shortening of this hard-core radius to 0.35a a p, 
but the speed of simulation is. Here e a p, and a a ^ are atomic interaction parameters, 
related to the well-depth, steepness and position, respectively, of the dispersion interaction. 
Second, a modified Axilrod- Teller term is added for all molecular triples. This is the triple- 
dipole dispersion correction to the van der Waals interaction which was first published by 
Axilrod and Teller.— In the derivation by Axilrod, 20 spherical polarizabilities are assumed, 
but for anisotropically polarizable molecules, such as CO2, his result does not hold. To 
investigate this case, we quickly recapitulate Axilrod's derivation^ - where we dispose of the 
assumption of spherical polarizabilities. 

For each of the three molecules I, m and n, we define a local coordinate system where 
the z-axis is normal to the plane of the molecular centers, the x-axis parallel to the bisector 
of the angle spanned by the two other molecules and the j/-axis mutually orthogonal to the 
x- and z-axes. Hence, the z-axes all coincide between the three local coordinate systems, 
but the x- and y-axes need not. Assume now that the electronic structure of each molecule 



is independent and described by a wavefunction that factorizes into separable x-, y- and 
^-contributions in its local coordinate system, i. e. for molecule I, 



TP l (x l ,y l ,z l ) = X{x l )Y{y l )Z{z l ) (9) 

We write the third-order perturbation correction to the groundstate energy, Wq", which in 
Axilrod's notation is (Eq. [5a] in Ref. Q), 

TJI TJI TJI 

w'" V V 0j jk k0 do) 



, in r hn v-y-i t ~\ f vi /-\ I / \ ~v v i /-\-n-f- /" 



Here H'- k is the matrix element of the dipole perturbation for states j and k, and W® is the 
energy of state j. Invoking the closure approximation,— this equation may be approximately 

recast as 

W °" H'^H'jA = «dtop,8(i, m, n) (11) 

j k 
j^kfi k^O 



where 



v' = ((Wf-W °)(W° k -W°)} jk (12) 



and (. . .)jk signifies the arithmetic average over j and k. Eq. (ITT]) serves to define the 
three-body correction to the dispersion energy that we will use. 

Formally, the set of matrix elements {Hj k } covers all possible excited states but we shall 
assume contributions to the sum only from the first excited orbital of each symmetry for each 
molecule. That is, the three lowest excited states of the arbitrary molecule I are assumed 
to be X*(xi)Y(yi)Z(zi), X(xi)Y*(y l )Z(z l ) and X(xi)Y(yi)Z*(zi) where the asterisk denotes 
the next higher-energy orbital. Because they share a common orthogonal z-axis, only mixed 
excited states for x- and y- components between the molecules contribute to the sum over 



states. Hence, the matrix e 
given in Table I of Ref. 
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ements for the sequences of possible excitations are exhaustively 
and the explicit form of these matrix elements is provided in 
Eq. (29) of the same reference, except that the common factor M 2 is no longer applicable 
because the transition diople moment is no longer the same for the different components. 
Instead, given the two molecules / and m, excited in their x and y orbitals, respectively, we 



write the corresponding matrix element (cf. Eq. [29] in Ref. |20|) 



{x m ,yi) = (M Em M w /(47rei2^))(2cos^7 m sin^ - sin -7™ cos ^) (13) 



where M Xm denotes the expectation value of transition dipole moment along the x-axis of 
molecule m and 7 m the angle defined by the molecules I, m and n with its the apex in molecule 



m. All the other matrix elements follow by analogy with Ref. |20|. We have written Eq. ( II 3p 
in a general form with e being the permittivity of the medium in which the molecules are 
dispersed. It is most reasonable to take this as the permittivity of free space. Since because 
of the very short-range nature of the interaction (it tapers off as the inverse ninth power 
of distance) it is unreasonable to assume that a homogeneous medium of CO2 molecules 
can be accommodated between the interacting molecules. The final approximation is to 
replace the M-factors by the square-root of the corresponding polarizabilities. Collecting 
the constants of proportionality in the common prefactor of Eq. (ITTj) . which is then seen to 
have dimensions of reciprocal energy, we treat it as a fitting parameter and denote it by 1/v 
proper. 

With anisotropic polarizabilities, the sum over states in Eq. (fTTj) does not simplify to the 
simple form given in the original references' 1 ^ and the complicated closed-form expression 
will not be reproduced here. In any case, since it involves sines and cosines it is not optimal 
from a computational point of view; it is much more efficient in terms of the total number 
of floating-point operations to calculate the truncated sum over states directly. To this end, 
the half-angle formulas are used to rewrite the matrix elements, such as the one in Eq. (fT3|) . 
in terms of dot products and square-roots, which are much more efficient in terms of CPU 
cycles than trigonometric functions. The polarizabilities are calculated, like before, as the 
interpolation 

a(0) = a± + I cos(6 l )|(a|| - a±) (14) 

where 9 is, once again, the angle to the molecular axis. The limiting polarizabilities are 
taken to be the same as the static ones. 

In total, after self-consistent solution of the induced dipoles, the energy of interaction 
among N molecules is given by 

N N N N N N 

1=1 m>l a£l (3£rn 1=1 1=1 m>l n>m 

The analytical gradient of this expression is very involved, with the chain-rule giving factors 
proportional to the gradient of the polarizability. Consequently, when the gradient has been 
needed, for instance, in energy minimization, it has been calculated numerically using the 
finite-difference approximation. 



C. Parametrization strategy 



1. Gas-phase properties 

A number of parameters have not been optimized, but simply assigned from plausible 
experimental values in the literature. Thus, the bond length is fixed at 1.161 A, midway 
between published values of 1.160 A and 1.162 A by experimental groups,— 1 ^ and the partial 
charges on the atoms are chosen to reproduce the experimental quadrupole moment.— & 
To further reduce the number of free parameters, the charge width of the oxygen atom, 
ro, was set equal to 0.610 A, the value of the oxygen atom in GCPM water ,— and, the 
corresponding quantity for carbon was scaled according to the ratio of the a parameters 
(vide infra). These parameters were not optimized. Moreover, I introduced the additional 
constraint of = 7^ = 7 Q /3 and for the remaining parameters, the following "mixing 
rules" were adopted for unlike interactions of the modified Buckingham exp-6 potential, 

e a/ 3 - - — — — (16) 

= ~ (&aa + CTpp) (17) 

Eq. (I16p can be justified with appeal to the London formula,-^ in which the harmonic average 
of the ionization energy is taken. Normally, it is the geometric average of the polarizabilities 
in said equation that lends its mathematical form to the combining rule for the e-parameters. 
However, the precise form of the mixing rules are of a secondary concern, as they serve 
mainly to reduce the parameter space that has to be fit, and provided the atomic interaction 
parameters do not turn out to be very different from each other, the result will be insensitive 
to reasonable choices of mixing rules. 

A test set of potentials with predefined a values, covering a broad range, but with the 
ratio, 0"c/o"o of the carbon sigma value, ctq, and the oxygen sigma value, <To, conserved 
at 1.0483 were investigated. This value was arbitrarily chosen early in the development of 
the model and never subjected to revision. For each such class of IPs, manual tuning of 
the e-parameters in trial-and-error fashion was made to obtain a good fit for B%(T) against 
experimental data and a reasonable binding energy and geometry of the dimer structure. 
B 2 (T) was calculated from the Mayer sampling method^ over 5 x 10 8 Monte Carlo steps. 
Typically, with this number of Monte Carlo steps, the resulting uncertainty, calculated from 
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the block average method— with blocks of 10 5 cycles, is of the order of one percent. The 
induced dipoles were iteratively solved to within a tolerance of 3 x 10~ 10 D at all times. In 
the course of this parametrization, it was found that the 7-parameter plays a most crucial 
role. This parameter controls the steepness of the modified Buckingham potential, without 
affecting well depth or location. Small values lead to a flat minimum and B 2 (T) turned out 
to be very sensitive to this parameter. Thus, it was found that, to fit B 2 (T), this parameter 
had to be increased from its initial estimate of 12.75 (taken from GCPM water—) to 15.50. 

Later in the development of the model, it was found necessary, with respect to the liquid 
densities, to include also the three-body dispersion in the energy expression. Because the 
^-parameter is completely independent of all dimer properties, including B 2 (T), it was fitted 
independently so that B%(T) coincided with the data of Dushek et al.,— It was not possible to 
reproduce, with this single parameter, also the data of Holste et al.,— something which lends 
credence to the former measurements. The value obtained, v = 2.82 x 10 4 K, is reasonable 
in that it is of the same order of magnitude as the Axilrod- Teller coefficient obtained for 
argon, for which this coefficient is 2.1 x 10 4 K in the appropriate units.— 



2. Bulk simulations 



The fluid densities were extracted from isothermal-isobaric Metropolis Monte Carlo 
simulations^ for an ensemble of 200 molecules in periodic boundary conditions of cubic 
symmetry over M = 2.5 x 10 7 steps, run in parallel over five independent Markov chains. 
Standard errors of the mean were estimated from the block average method^ with \/M 
blocks. 

For the modified Buckingham potential, with neglect of interactions beyond the cutoff 
the energy was corrected by 

f^LRC = tlgfe + 4n^ c + ^lrc (18) 

where mlrc — f Jr°° 4 7rr ' 2 f4xp-6( r )dr and superscripts indicate between which two atom 
types the interaction is computed. Here p is the number density of molecules. The integral 
in question can be analytically computed, which yields 



2nep 

mlrc = : 



6/7 



— (r c7 +2r c7 a + 2a ) exp (7 - — j - — 



(19) 
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For the Axilrod- Teller terms, in light of their very strong distance dependence, no long-range 
correction was deemed necessary. 

For the electrostatic interaction, a long-range correction was introduced for the induced 
dipole using the Onsager expression for the dipole reaction field^ with the dielectric con- 
stant taken from the Clausius-Mosotti equation fitted against experimental measurements 
of the dielectric constant.— This may seem to be very approximate but was done for two 
principal reasons. First, lattice summation techniques such as the Ewald sum^ introduce 
assumptions of periodicity which are unsuitable for the uncorrelated nature of the molecules 
in a disordered phase. Second, the quadrupolar interaction terms decay sufficiently fast 
with distance so as to make their sum absolutely convergent. On the assumption of un- 
correlated molecules beyond the cutoff radius, the long-range correction vanishes. Last but 
not least, the exact same procedure was used by Paricaud et al. in their simulations of 
GCPM water.- 14 As a technical side-note, it should be pointed out that Ewald summation 
for Gaussian charges is complicated by the fact that the required Fourier transform for the 
reciprocal space sum is not analytically tractable. Using the reciprocal space sum for point 
charges, with due corrections in direct space and a long cutoff radius, is an alternative but 
inefficient. 

Each trial move consisted of randomly displacing and rotating from one up to four 
molecules. Interaction cutoffs were introduced at half the box length. For the three-body 
dispersion, this was interpreted to mean that all interacting molecules had to be within 
cutoff of each other. Simple mixing^ with a mixing factor of 10% was used to enforce 
and speed convergence. Because of increased computational load, the induced dipoles were 
solved to within a tolerance of 3.4 x 10~ 5 D per molecule. This is still less tolerant than 
many other reported simulations of polarizable IPs: Ren and Ponder— report simulations 
on liquid water with their AMOEBA force field where the dipoles are solved to within 10 -2 
D precision; Paricaud et al.. 14 solve the induced dipoles to within 5 x 10~ 5 D in their sim- 
ulations of GCPM water. However, further relaxation of the tolerance is unnecessary as 
the greater part of the simulation time is not spent on the self-consistent solution of the 
elctrostatic field equations, but (about 70%) on the evaluation of J2imn M disp,3(^, m, n). This 
computation can be significantly sped up with shorter interaction cutoffs. 

Two types of bulk simulations were performed. A series of iVpT-simulations to determine 
the density and constant-pressure heat capacity of the model fluid and iV^/T-simulations to 
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determine the radial distribution functions of the model fluid. Together with the calculation 
of the virial coefficients, these served to parametrize the model. All parameters except 
for v were decided for the potential model lacking the % term, which has no effect on 
either the dimer binding energy, its geometry or the second virial coefficient. Selection 
among candidates of this pairwise dispersion model was effected by comparing densities 
from the iVpT-simulation with experiment, although the long-range correction was not fully 
developed at this time. Subsequent introduction of the long-range correction lead to an 
increase of the computed densities, later corrected by the introduction of the three-body 
dispersion interaction. The final parameters, extracted from these tests, are listed in Table 
i 



III. RESULTS AND DISCUSSION 
A. Virial coefficients 

Not only have calculations for the parametrization of the GCPCDO IP been carried out, 
but in the investigations leading up to its formation, I also investigated some other COo 
IPs. The VLE envelopes of these models have all been thoroughly investigated in Ref. |6, 
in which it was established that the ZD potential exhibits near-perfect experimental agree- 
ment in this property. While this excellent agreement might not be reproducible in every 
aspect,— the critical temperature and density cannot lie very far from their experimental 
counterparts nonetheless, because the relative deviations between the coexistence densities 
reported in Ref. O and Ref. [l2 decrease and approach the experimental values when nearing 
the experimental critical temperature. Moreover, the variation in B2(T) is not very great 
between the models; except possibly for the TraPPE IP— which — it is interesting to point 
out — is that of the non-polarizable IPs of MSM-type, which exhibits the best experimental 
agreement for both B 2 (T) and B 3 (T). For the GCPCDO IP, having had this as one of its 
goals in the parametrization, the fit is very good with the relative error never exceeding 5%. 
The results are reported for a select number of temperatures in Table [Til 

It does seem surprising that IPs that fare very well in reproducing VLE properties should 
fail so remarkably at reproducing the much "simpler" property B 2 (T). If B 2 [T) is overesti- 
mated, then a compensating underestimation of B 3 (T) seems very likely. Direct calculation 
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TABLE I. Interaction parameters and experimental observables for the GCPCDO model and, 
where appropriate, experimental reference values. 



Parameter 



GCPCDO 



Exp. 



do / A 
a c I A 
eo /K 
e c /K 
v I K 
7o / A 
7c / A 
to I A 

r c / A 

qo I e 
<?c / e 
a ± / A 3 
ay / A 3 
Bond length / A 

Quadrupole / (e A 2 ) 



3.347 
3.193 
67.72 
71.34 
2.82 x 10 4 
15.50 
15.50 
0.6100 
0.5819 
-0.3321 
0.6642 
1.95 
4.05 
1.161 

-0.90 



1.929* 
4.038.1 
1.160^ 
1.162°. 
-0.85f. 
-0.90^ 
-0.96 c 



a Ref. 
b Ref. 
c Ref. 
d Ref. 
e Ref. 
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25 



of B 3 (T) for the computer models seem to confirm this. However, contrary to the case of 
B2(T), good quality experimental data are very hard to find for B^{T). Different authors 
report widely different results, over the same temperature range. In the narrow range around 
the critical temperature, however, both Holste et al— and Dushek et al— report measure- 
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TABLE II. Second virial coefficients for the ZD, MSM, EPM2, TraPPE and GCPCDO IPs, as 
well as experimental results, at a select number of temperatures. For the computed results, brack- 
eted numbers indicate the estimated standard error of the mean in the last digit from the Mayer 
sampling 27 Monte Carlo integration. 





B2 1 (cm 3 


mol 1 ) 










T / K 


ZD 


MSM 


EPM2 


TraPPE 


GCPCDO 


Exp. 


220.00 


-200.6(6) 


-204.7(6) 


-206.6(7) 


-222.0(7) 


-247.2(9) 


-247.501 
-247.521 


240.00 


-168.1(5) 


-171.5(6) 


-170.9(6) 


-183.5(6) 


-202.8(9) 


-202.831 
-202.131 


260.00 


-141.7(5) 


-145.0(5) 


-144.7(5) 


-153.7(6) 


-169.8(9) 


-168.921 

— IDO.Z (_ 


280.00 


-120.6(5) 


-123.3(5) 


-122.1(5) 


-129.3(5) 


-143.7(7) 


-142.701 
-142.il! 


300.00 


-104.4(4) 


-106.2(4) 


-105.0(4) 


-110.3(5) 


-123.5(6) 


-121.701 
-121.351 


340.00 


-78.0(4) 


-79.5(4) 


-78.0(4) 


-81.9(4) 


-91.3(4) 


-90.571 


423.15 


-43.5(3) 


-44.8(3) 


-43.8(3) 


-45.6(4) 


-52.2(4) 


-51.251 


448.15 


-36.6(3) 


-37.5(3) 


-36.5(3) 


-37.7(3) 


-43.9(3) 


-43.511 



a Ref. 
b Ref. 
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ments which are in at least slight mutual concordance. In Table ITTTj. these measurements are 
reported and compared with predictions from the IPs. As is evident, the pairwise additive 
IPs underestimate B 3 (T) across the whole temperature range. 



B. Volumetric properties 

To investigate the properties of the many-body potential with more than just three bod- 
ies, the density and the heat capacity at constant pressure, both readily extracted from the 
NpT simulations, serve as indicators. These results are summarized in Table IIVI with ex- 
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TABLE III. Third virial coefficients for the ZD, MSM, EPM2, TraPPE and GCPCDO IPs, as well 
as experimental results, at a select number of temperatures. For the computed results, bracketed 
numbers indicate the estimated standard error of the mean from the Mayer sampling 2 ? Monte 
Carlo integration. 



T / K 


I (cm 6 
ZD 


mol 2 ) 

MSM 


EPM2 


TraPPE 


GCPCDO 


Exp. 


280.00 


2920(30) 


2940(30) 


2910(30) 


3080(40) 


5140(60) 


56361 














5165^ 


300.00 


2820(20) 


2860(20) 


2922(20) 


3060(30) 


4790(60) 


4927f. 














4753^ 


320.00 


2690(20) 


2760(20) 


2700(20) 


2870(20) 


4460(50) 


44231 














4360^ 


340.00 


2530(20) 


2570(20) 


2560(20) 


2680(20) 


4046(50) 


3996^ 



a Ref. 
b Ref. 
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perimental data. It is surprising that the agreement with experiment is considerably better 
at the higher densities. This increasing discrepancy between theory and experiment can 
be tentatively attributed to the fourth virial coefficient, B±(T). It is clear that in order to 
explain these results, B±(T) (or possibly higher virial coefficients) must rise quicker with 
temperature for the GCPCDO model than for its experimental counterpart. Unfortunately, 
quality experimental values of B±(T) are not known but it is worth pointing out that if one 
applies the virial equation of state truncated after B%(T), the computed densities at 290 K 
and 310 K turn out to be 0.915 ± 0.018 g / cm 3 and 0.791 ± 0.015 g / cm 3 , respectively, 
at 200 bar pressure. These values are indeed very close to the simulated values reported 
in Table HVl and indicate that B±(T) is overestimated and close to zero. However, because 
at the higher pressure, the truncated virial equation of state predicts densities higher than 
the simulated ones, it is clear that higher-order coefficients must be compensating for errors 
in B±(T). This is not surprising, at high densities, the steric effects become the dominant 
mode of interaction. Sadly, the precise causes of these deviations in the virial coefficients 
is unknown. Unlike the other models investigated in this work, we have come some way in 
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correcting B2(T) and B$(T) to their correct experimental values. However, we may safely 
say that the interactions of the CO2 molecule are not as simple as one may at first suppose. 

Also shown in Table HVl is the constant-pressure heat capacity which was calculated from 
the fluctuation formula,-^ 

C P = j^r 2 ((H 2 ) - {Hf) (20) 

where H = U + pV + 5kT/2 is the enthalpy, the last term being the classical kinetic 
contribution of a linear rigid body and k the Boltzmann constant, p the pressure, U the 
potential energy, V the volume, N the number of molecules and T the temperature. For 
a completely fair comparison with the experimental values, also the internal vibrational 
degrees of freedom should be included. Assuming harmonic behavior, for each normal mode 
of frequency v this quantized harmonic contribution to the heat capacity is then 

^KS'MS -1 ) (2i> 

where h is the Planck constant. Taking into account the experimental frequencies^ 4 ^ of the 
four harmonic normal modes of the CO2 molecule, this term is added to C p and reported as 
the corrected values in Table IIVI Not surprisingly, this expression compares favorably with 
the experimental C p extrapolated to vanishing density— For instance, the experimentally 
ascertained intramolecular contribution is 5.74 J / (K mol) at 250 K, whereas from Eq. ( f2~TT) 
one has 5.78 J / (K mol). At 310 K, the experiments indicate 8.58 J / (K mol) and 
from Eq. (pTI) we have 8.68 J / (K mol). It is computationally too demanding, at present, 
to include the quantized vibrations in the bulk simulation of the GCPCDO IP and the 
approximation of separable internal and external degrees of freedom is expected to be fair. 

The general overestimation of the heat capacity, even before the correction for intramolec- 
ular degrees of freedom, is due, at least in part, to the assumption of classical translational 
and rotational degrees of freedom. Especially at high density, free rotation and translation 
is not possible and the rotational and translational degrees of freedom are in effect partly 
quantized librational modes. Unlike the intramolecular degrees of freedom, these are highly 
coupled and there exists no viable computational approximation for their contribution to 
the heat capacity. The very large overestimation of the heat capacity at 310 K and 200 bar 
cannot, however, be attributed to this effect alone. Moreover, at this state point the dis- 
crepancy in density between real CO2 and GCPCDO is so large that a more fair comparison 
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TABLE IV. Calculated and experimental fluid properties at 200 and 800 bar pressure and tem- 
peratures of 250 to 310 K. Heat capacities assume classical contribution of 3k for the rotational 
and translational degrees of freedom. The corrected heat capacities include the heat capacity of 
the quantized internal degrees of freedom within the harmonic approximation. For the calculated 
results, numbers in parentheses indicate the estimated standard error of the mean in the last digit. 



p I bar 



T / K 



P I (g / cm 3 ) 



C p / (J / [K mol]) 







Calc. 


Exp.^ 


Calc. 


Corr. 


Exp.^ 


200 


250.0 


1.100(4) 


1.105 


96(4) 


102(4) 


83.3 


200 


270.0 


1.010(5) 


1.032 


87(3) 


95(3) 


86.0 


200 


290.0 


0.902(7) 


0.951 


91(3) 


99(3) 


90.4 


200 


310.0 


0.79(1) 


0.856 


109(5) 


118(4) 


97.7 


800 


250.0 


1.216(2) 


1.211 


76(3) 


82(2) 


74.6 


800 


270.0 


1.159(2) 


1.165 


71(2) 


78(2) 


73.7 


800 


290.0 


1.107(2) 


1.118 


72(2) 


80(2) 


73.0 


800 


310.0 


1.057(2) 


1.073 


68(2) 


77(2) 


72.4 



Ref. 
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(as relates to C p ) is with the experimental C p at 0.79 g / cm 3 density, which is 42 119 J / (K 
mol), not too far off from the computed value. 

No IP for the fluid can be deemed satisfactory if unable to reproduce the experimental 
fluid structure. Accordingly, at the density and temperature of the neutron-diffraction 
experiments of Cipriani et al.,— the atomic pair distribution functions (PDF), g(r), have been 
computed, and these are shown in Figures [1] and |2] for two different thermodynamic states. 
What is experimentally ascertained, however, is not the individual, atomically resolved g(r), 
but the superimposed effect from all atomic scatterers. Taking account not only of the four 
times greater abundance of C-0 and 0-0 vectors than of C-C vectors between molecules, 
but also of the different propensity toward neutron scattering, the following formula was 
used to calculate the superimposed PDF,— 

G(r) = 0.403(?oo(r) + 0A6Ag co (r) + 0.133£ C c(r) (22) 

In terms of the position and number of peaks, these G(r) are seen to be in reasonable-to-good 
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FIG. 1. The individual atomic PDFs, g(r), of Eq. (|22p for carbon-carbon (long-dashed line), 
carbon-oxygen (short-dashed line) and oxygen-oxygen (dotted line) at 312 K and 0.83 g / cm 3 for 
the GCPCDO IP as well as their superposition (full line), weighted by occurrence and scattering 



propensity, G{r). Squares are experimental results from Ref. 
contributions are shown for the computed results. 



43J. For clarity, no intramolecular 



agreement with the experimental results. 

At 240 K and 1.09 g / cm 3 , the carbon-oxygen distribution function clearly shows more 
structure at short range, than at 312 K and 0.83 g / cm 3 , where the lack of orientational 
correlation in the fluid is also apparent in how quickly the atomic carbon-oxygen and oxygen- 
oxygen PDFs decay to unity. Still, however, the carbon-carbon PDF exhibits a slight peak at 
around 7.5 A indicating a weak second coordination shell, albeit of random order in molecular 
orientation. The PDFs indicate that the fluid is slightly overstructured at the higher density, 
where the first peak is overestimated. The better agreement for the computed PDF at the 
low density is most likely due to the overall lesser contributions from many-body effects at 
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the first peak is overestimated at both 



this density. For the pair potential used in Ref. 
thermodynamics states. Last, it should be pointed out that in the fluid with flexible bonds, 
a general broadening of the peak structure is expected. This effect is at least responsible 
for the deviation seen at the very short distances, where the internal scattering vectors 
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FIG. 2. Same as Figured] but for 1.09 g / cm 3 and 240 K. 
contribute. In the rigid model, they are ^-functions and have been omitted for clarity. 

C. Clusters 

Having established the weakness of the GCPCDO model primarily in its B^T), i. e. four- 
body interaction, we turn to properties of the IP for which only three bodies contribute. 
Clusters like these offer excellent tests of the model, due to the availability of quantum- 
mechanical reference calculations. It also allows us to pinpoint more clearly the role of 
many-body effects in the interaction potential. 

1. Dimer and trimers 

Both experiments^ and ab initio simulation^- 3 - agree that the equilibrium dimer struc- 
ture is of C*2h symmetry, with the ab initio simulations indicating that there is a saddle point 
of Civ symmetry. The GCPCDO IP reproduces these two dimer states very well, as shown 
in Table |V] As for the geometry of the dimer configurations (see Figure |3j) , it is neither 
better nor worse than the simpler ZD IP,- but when it comes to the binding energy of the 
two states, it is markedly superior when judged against the ab initio IPs: the binding energy 
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FIG. 3. Schematic illustration of the dimer with definitions of angles and distances. Each molecule 
is represented by a stick. All atoms lie in the same plane. 

is underestimated by less than 4% for the minimum and overestimated by less than 1% for 
the saddle-point, whereas the ZD IP underestimates the binding energy in both cases by 
about 15-20%. This gross underestimation of the dimer binding energy helps explain the 
overestimation of B 2 (T) for the ZD model (vide supra) but is nevertheless surprising be- 
cause for non-polarizable potentials the general trend is an overestimation of dimer binding 
energies. 

Because the GCPCDO model is so successful at capturing the dimer binding energy, it is 
interesting to test it across a broader range of conformations. Hence, I show in Figure H] the 
potential energy function of this model compared to the accurate data of the angular fit of 
the symmetry-adapted perturbation theory^, calculations of Bukowski et al. ,— the BUK IP, 
for radial displacements. Both the energy minimal pathway separating the two molecules 
and the path along the C 2v transition state are shown. For the energy optimal pathway, 
the optimized angles of the two molecules at each separation are shown in Table IVIl The 
two molecules, for both the GCPCDO and BUK IPs, prefer a slipped-parallel conformation 
at close range, but eventually prefer the T-shaped geometry of the minimum at long range. 
The transition is noticeable as a slight trough in the dissociation curve at around 4 A and 

19 



TABLE V. Equilibrium geometry and well-depth energy of the CO2 dimer for the GCPCDO 
model, the ZD model, and BUK, the angular ab initio potential energy surface of Bukowski et al£ 

U refers to the potential well-depth at the specific conformation. See Figure [3] for the definitions 
of the geometric quantities. 

GCPCDO ZD BUK? Exp.^ 

Minimum (C2/J 

U I K -675.4 -548.7 -696.8 

01,02 /deg 55 - 5 56 - 5 59 -° 57 - 96 

R I A 3.64 3.64 3.54 3.60 

Saddle-point {Civ) 

U I K -596.6 -508.7 -593.2 

01 / deg 90.0 90.0 90.0 

2 / deg 0.0 0.0 0.0 
R I A 4.16 4.17 4.14 



a Ref. 
b Ref. 



2 

4G 



is quicker for the BUK IP with an earlier onset. 

Another interesting property of the model, which cannot be answered by the BUK IP, 
is the total dipole of the Ci v configuration. Because the electric field gradients are very 
inhomogeneous close to the molecule, it might be suspected that only allowing the center 
atom to polarize is artificially deflating the induced dipole, and like this introducing er- 
rors in the short-range interaction. For GCPCDO, the dipole is predicted to be 0.19 D 
at the transition-state separation. Running calculations with Gaussian 03 ,— the suspi- 
cion is confirmed as calculations of the dipole moment in the identical configuration at the 
CISD/aug-cc-pVDZ level of theory, predict a dipole moment of 0.21 D. The dipole moment 
is hence underestimated by 10% in the transition state configuration. 

Because of its many-body nature, it is interesting to test the GCPCDO IP on the simplest 
cluster for which many-body effects contribute, i. e. the trimer. Consequently, energy 
minimized structures of the trimer have been located. The two most stable conformations 
are shown in Figures [5] and H2 the specific data on each are summarized in Table IVHI 
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FIG. 4. Dependence of energy of interaction, U(R), on separation, R, for two different paths of 
approach. Squares are for the energy minimal path of the GCPCDO IP and the full line is the 
analogous result for the BUK IP. Circles are for the Civ conformation of the GCPCDO dimer, and 
the dashed line is for the BUK result. 

Both of these trimer conformations have been observed spectroscopically, with the planar 
C*3/j conformation being slightly more abundant.— 1 ^ In terms of relative energies, however, 
neither of GCPCDO or BUK predict the right ordering, but there are two general remarks to 
be made. The first one is that the inclusion of many-body effects, clearly levels the difference 
between the two states, the difference in the GCPCDO prediction being less than 0.1 K.— 
That many-body effects would alleviate the problem was hinted at already by Bukowski and 
coworkers^ in their discussion of this problem and they argued using single-point calculation 
from higher-order quantum chemical theory that this was the case. Second, the possible 
role of the zero-point vibrational quanta has to be kept in mind. A first-order estimate of 
this effect is provided by the harmonic zero-point energy. Indeed, as indicated by Bukowski 
and coworkers,- inclusion of this energy for the BUK IP decreases the difference between 
the states. Carrying out the same analysis for the GCPCDO IP, however, we find that 
theory is brought into qualitative agreement with observation. As discussed by Bukowski 
and collaborators,- the harmonic approximation is very strained in the CO2 trimer but 
correct evaluation of the zero-point vibrational energy necessitates a numerical solution of 
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TABLE VI. The optimum values of the angles 0\ and 02 as a function of the separation R be- 
tween two molecules as predicted by the GCPCDO IP and the BUK IP. The angles are defined 
geometrically in Figure [3j 





GCPCDO 




BUK 




R /A 


0i / deg 


6o I deer 


0i 1 deg 


02 / deg 


3.0 


71.7 


71.7 


66.4 


66.4 


3.3 


62.3 


62.3 


62.2 


62.2 


3.5 


58.0 


58.0 


59.5 


59.5 


3.7 


54.5 


54.5 


57.0 


57.0 


3.8 


53.0 


53.0 


46.1 


65.2 


3.9 


42.1 


61.3 


37.1 


71.1 


4.0 


30.8 


70.2 


28.6 


75.9 


4.1 


19.6 


77.9 


18.9 


80.9 


4.2 


3.7 


88.0 


0.0 


90.0 


4.5 


0.0 


90.0 


0.0 


90.0 



the Schrodinger equation. Future code development will remedy this situation. Even so, I 
believe that these results are indicative of the qualities of the GCPCDO IP. 

2. Water complexes 

Because of its transparent physical form, the GCPCDO IP can, through the adoption 
of suitable "combining rules", be interfaced with other IPs. As a first test of the feasibil- 
ity of this approach, I have calculated the binding energy and molecular geometry of the 
[H 2 0— C0 2 ] complex using the successful GCPM water— for the water moiety. In addition 
to the combining rules of Eqs (|T6|) and (ITT)) , the 7-parameters were calculated as 

lxy=^{lx + ly) (23) 

For comparison purposes, the complexation of ZD CO2 and TIP3P water™ serves as an 
indicator of the effect of neglected polarization. These results are summarized in Table 
IVIIIl It is important to point out that the potential energy surface of this complex is very 
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FIG. 5. Schematic illustration of the C^h trimer. All atoms lie in the same plane. Molecules A, B 
and C are all identical by symmetry. 




View from the top 



c 

j View from the side 

\P 




FIG. 6. Schematic illustration of the C2 trimer. Molecules A and B are identical by symmetry and 
the distance between their centers is R. The perpendicular distance from the center of the line 
joining their centers to the center of molecule C is P. 
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TABLE VII. Energy well depth, U, equilibrium geometry and harmonic zero-point vibrational 
energy, Eq, for the two most stable trimers predicted by the GCPCDO IP and the BUK IP. 
Because of uncertainties in the numerical Hessian, Eq is rounded for the GCPCDO. The geometric 
variables are defined in Figure [6] for the C2 minimum and in Figure [5] for the C^h minimum. 

GCPCDO BUKf. 

C2 minimum 

U I K 
E /K 
R/k 

Pj A 

a / deg 
/3/deg 
7 / deg 

C^h minimum 

U I K 
E /K 
R/k 
P I deg 

a Ref. 2 

shallow near the minimum so precision is difficult, and two different ab initio minimum 
energy geometries have been reported in the literature. Danten and collaborators^ find 
that the energy minimum has C s symmetry from counterpoise-corrected MP2 calculations 
with the aug-cc-pVTZ basis set. However, this is contested by both experimental^ 1 ^ and 
more high-level ab initio counterpoise-corrected calculations at the CCSD(T) level of theory 
and the same basis set which predict symmetry of the complex.— It is interesting to note 
that the MP2 level of theory often overestimates correlation effects, consistent with the fact 
that simple pair potentials, such as the ZD& and TIP3P 52 IPs under the Lorentz-Berthelot 
mixing rules, predict a C<i v minimum for the complex. The intermixing of the GCPCDO 
IP with GCPM water 14 produces two, equivalent, global minima of C s symmetry, but the 
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-1849.4 
359 
3.73 
2.96 
12.1 
53.7 
152.8 

-1849.4 
318 
4.07 
40.9 



-1889.8 
343.6 
3.60 
2.89 
10.9 
53.4 
157.6 

-1819.6 
304.2 
4.04 
39.3 



difference in energy between them and the transition state of C% v symmetry connecting them 
is only about 7 K, i. e. less than 0.1% of the total interaction energy and should not be 
taken as great drawback of the model. 

It is a greater error that the total binding energy is underestimated by about 25%, a fault 
shared by the ZD/TIP3P combination as well. This is surprising, because in general non- 
polarizable IPs parametrized against bulk properties tend to overestimate cluster binding en- 
ergies. The explanation can partly be found in the rigid geometries of the GCPCDO/GCPM 
and ZD/TIP3P moieties. The results by Danten and coworkers^ indicate that the CO2 
molecule is bent slightly upon complexation, thus further inducing a dipole moment which 
increases the energy of interaction. However, the greater part of the explanation is to be 
found in the approximation that only the central atom polarizes. A QM calculation at the 
CISD/aug-cc-pVDZ level of theory for the C s minimum of the GCPM/GCPCDO complex 
indicate that while the total dipole moment of 2.19 D is in good agreement with the accurate 
value of 2.21 D, the individual components are poorly reproduced.— The QM calculation 
indicates that for the complex as a whole, the electric dipole moment along the C2 axis 
of the water molecule is 2.16 D, and that the perpendicular component is 0.465 D in this 
particular configuration. The GCPM/GCPCDO pairing, however, indicates 2.19 D and 0.14 
D, respectively, for these components. Further discrepancy is expected at the closer range 
indicated as the equilibrium distance by the ab initio calculations. In summary, as far as 
geometry and energy are concerned, the predictions of the GCPCDO/GCPM complex can 
hardly be considered an improvement over the simple empirical IPs, compared to ab initio 
calculations. 

Despite this small short-coming of the predictions for the [C0 2 — H 2 0] complex, it is still 
worthwhile to consider the [(C0 2 ) 2 — H 2 0] complex, also studied by Danten and coworkers,— 
because here the true many-body interactions start to play a role. Moreover, contrary to 
the case of the C0 2 trimer, where none of the moieties is dipolar, the water molecule carries 
a substantial dipole moment and the electronic induction effects are expected to be more 
pronounced. It must be pointed out that despite this being a three-body system, no Axilrod- 
Teller potential has been applied. The reason is that while GCPM water has a known 
polarizability, it has no Axilrod- Teller coefficient. Rather than impose one on the model, I 
have decided to judge it fairly according to its own merits, and these exclude a three-body 
dispersion interaction. The results are shown in Table IIXI The agreement is very good in 
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TABLE VIII. Minimum energy, U, and geometry of the [C0 2 — H 2 0] complex as predicted by the 
mixing of GCPCDO and GCPM polarizable IPs, the ZD and TIP3P non-polarizable IPs, or ab 
initio calculation by Danten et al. at the MP2/aug-cc-pVTZ level of theory, or by Garden et al. 
at the CCSD(T)/aug-cc-pVTZ level of theory. The binding energy of the transition state 
predicted by the GCPCDO model is reported for completeness. For the definition of the geometry, 
see Figure [3 





GCPCDO/GCPM 
C s minimum 


C2v tr. state 


ZD/TIP3P 


Danten^ 


Garden^ 


U 1 K 


-1036.4 


-1029.5 


-1048.1 


-1308.4 


-1358.7 


R/k 


3.06 


3.06 


2.93 


2.77 


2.81 


<P 1 deg 


17.0 


0.0 


0.0 


13.9 


0.0 


a 1 deg 


87.4 


0.0 


0.0 


88.0 


0.0 



a Ref. 
b Ref. 
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FIG. 7. Schematic illustration of the structure of the monohydrate complex. The C2 axis of the 
water molecule is marked. All atoms lie in the same plane. 
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TABLE IX. Characteristics of the global minima of the [(C0 2 ) 2 — H 2 0] complex predicted by the 
GCPCDO/GCPM IPs, the ZD/TIP3P IPs or the ab initio results of Danten et al. at the MP2/aug- 
cc-pVTZ level of theory. The oxygen of the water molecule is placed at the origin of a Cartesian 
coordinate system in which the y-axis coincides with the C 2 axis of the water molecule with the 
positive direction pointing toward the hydrogen atoms, the z-axis is normal to the molecular plane 
and the x-axis is orthogonal to the y— and z-axes. The centers of mass of the CO2 molecules are 
given with respect to this coordinate system, and the orientation of each molecule is expressed in 
the Euler angles a, f3 and 7 which denote counterclockwise rotation around the x-, y- and z-axes, 
respectively. U is the energy of that conformation. Indices 1 and 2 denote the two CO2 molecules. 
The work in Ref. 53 used flexible molecules. The coordinates reported here are rounded over these 
differences. 





GCPCDO/GCPM 


ZD/TIP3P 


Danten^ 


Point group 


c s 


c 2 


C s 


U 1 K 


-2738.5 


-5074.0 


-2918.7 


xx 1 A 


0.00 


-1.43 


0.00 


yi 1 A 


-2.61 


-2.11 


-2.52 


z x /k 


1.58 


-0.85 


1.11 


x 2 1 A 


0.00 


1.43 


0.00 


V2 1 A 


0.79 


-2.11 


0.42 


Z2/k 


3.80 


0.85 


3.80 


ax 1 deg 


53.7 


-62.8 


61.8 


Pi 1 deg 


90.0 


86.0 


90.0 


7i / deg 


90.0 


-82.3 


90.0 


a 2 1 deg 


-69.1 


62.8 


-55.0 


fa 1 deg 


90.0 


86.0 


90.0 


72 / deg 


90.0 


82.3 


90.0 



Rcf. 



53 



27 



terms of energy. Contrary to the [C0 2 — H 2 0] complex, for the [(C0 2 ) 2 — H 2 0] complex, there 
is no appreciable underestimation of the binding energy. This can probably be attributed 
to the larger fraction of C0 2 -C0 2 interactions in this complex. Moreover, both Danten et 
alM and I find that the minimum is of C s symmetry, but for the ZD^/TIP3P~ model, the 
global minimum is of C 2 symmetry in a conformation reminiscent of the C 2 trimer. A local 
minimum of C 2 symmetry is predicted by the GCPCDO/GCPM— model at about 208 K 
above the global minimum. Clearly, polarization changes the relative stability of these two 
conformations in favor of the C s one and this is captured both by the MP2 calculations of 
Danten et al— and by the present work. The peculiarities of the global minima predicted 
by the IPs are given in Table IIXI It is very clear that the many-body effects are responsible 
for the altered symmetry of the equilibrium structure. Also, the binding energy is vastly 
overestimated by the non-polarizable IP pair. 

IV. CONCLUDING REMARKS 

A new, polarizable IP for C0 2 has been introduced and shown to be in excellent agreement 
dimer properties and excellent-to-passable agreement for the bulk phase. Classical non- 
polarizable IPs that reproduce bulk phase properties well have been shown inadequate for 
the description of B 2 (T) and also B 3 (T). That many-body effects should not be ignored for 
the CO2 molecule is evidenced by the qualitative experimental agreement that is achieved 
for the stability of the two trimer conformations, something which not even the ab initio 
potential energy surface of Bukowski et al- manages. Further corroboration of the model 
is provided by calculations on the water complexes, where especially the [(C0 2 ) 2 — H 2 0] 
complex is in good agreement with ab initio calculation at the MP2/cc-aug-pVTZ level of 
theory— Absence of polarization changes the symmetry of this complex. 

A tough test for all of the molecular C0 2 potentials investigated is the prediction of 
the virial coefficients. The GCPCDO model handles B 2 (T) and B 3 (T), but only because 
of design, and fails at B^(T) and up. This can still be considered an improvement over 
the classical, non-polarizable models for which most probably none of the virial coefficients 
beyond the ideal gas term are in agreement with experiment. Clearly, the interaction among 
C0 2 molecules is more complicated than a simple pairwise sum over atomic charges and 
Lennard- Jones terms. However, it is also more complicated than self-consistent solution of 
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induced dipoles and triple-dipole dispersion interaction. Nevertheless, for systems of three 
bodies or less, it seems to be highly satisfactory, meaning that the remaining errors relate to 
many-body effects beyond the third. On the precise causes of the remaining errors in the IP 
can only be speculated and future computer experiments may provide the answer to what 
the mechanisms are. 

It is in this light that it must be kept in mind that even if in comparison with experiment, 
the GCPCDO model, with a few exceptions, rests more on qualitative concordance than on 
many digits of precision, this is the first many-body molecular IP for the CO2 molecule to 
be developed and that many areas of inquiry remain to explore. One obvious and immediate 
improvement to the model would be to distribute the induced dipole over all atoms for a 
better short-range description of the induced electrostatics. 

The Fortran 90 source code for the GCPCDO energy subroutines is available upon re- 
quest. 
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